In-plane force fields and elastic properties of graphene 
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Bond stretching and angle bending force fields, appropriate to describe in-plane motion of 
graphene sheets, are derived using first principles' methods. The obtained force fields are fitted 
by analytical anharmonic energy potential functions, providing efficient means of calculations in 
molecular mechanics simulations. Numerical results regarding the mechanical behavior of graphene 
monolayers under various loads, like uniaxial tension, hydrostatic tension, and shear stress, are 
presented, using both molecular dynamics simulations and first principles' methods. Stress-strain 
curves and elastic constants, such as, Young modulus, Poisson ratio, bulk modulus, and shear mod- 
ulus, are calculated. Our results are compared with corresponding theoretical calculations as well 
as with available experimental estimates. Finally, the effect of the anharmonicity of the extracted 
potentials on the mechanical properties of graphene are discussed. 



I. INTRODUCTION 

Graphene consists of a two-dimensional (2D) sheet of 
covalently bonded carbon and forms the basis of both 
ID carbon nanotubes, 3D graphite but also of impor- 
tant commercial products, such as, polycrystalline car- 
bon (graphite) fibres. As a defect-free material, graphene 
is predicted to have an intrinsic tensile strength higher 
than any other known material [l[ and tensile stiffness 
similar to values measured for graphite. Graphene over 
the years and even prior to its isolation has been an ideal 
material to model as far as its mechanical properties are 
concerned. In particular, the elastic moduli of single layer 
graphene and its elastic response, have been a subject of 
intensive theoretical research in recent years and quite 
different approaches have been employed 0413 ■ For ex- 
ample, several groups have performed first principles' cal- 
culations @ .other used empirical potentials for atomistic 
simulations [4, 8, 9], and also tight-binding methods have 
been employed @, Q . As it is evident from the literature 
survey there is a large discrepancy of values regarding 
the stiffness of monolayer graphene and values ranging 
from 0.5 TPa to 4 TPa have been proposed depending 
on the methodology pursued in each case. 

From the experimental point of view, recent exper- 
iments have indeed confirmed the extreme stiffness of 
graphene of 1 TPa and provided an indication of the 
breaking strength of graphene of 42 N/m (or 130 GPa 
considering the thickness of graphene as 0.335 nm) [i~3| . 
These experiments involved the simple bending of a tiny 
flake by an indenter on an AFM set-up and the force- 
displacement response was approximated by considering 
graphene as a clamped circular membrane made by an 
isotropic material. An alternative way to assess how 
effective graphene is in the uptake of applied stress or 
strain (uniaxial or biaxial) is to probe the vibrational 
characteristics of certain interatomic bonds upon load- 
ing. In particular, Raman spectroscopy has been proved 
very successful in monitoring the mechanical response of 
certain Raman active modes upon the application of ex- 



ternal mechanical stress/strain. In axial tension, a linear 
relationships between Raman frequency and strain was 
established regardless of the geometry of the monolayer 
graphene up to strains of about 1.5 % [T3 - [T7j . 

In this work, we present appropriate empirical force 
fields, derived from first principles' calculations, to de- 
scribe bond stretching and angle bending interactions in 
graphene. Analytical nonlinear potentials are provided 
for these interactions, in order to efficiently implemented 
in atomistic molecular dynamics or Monte Carlo simula- 
tions. A big advantage of atomistic models compared to 
the more accurate first principles' methods is their com- 
putational efficiency which makes them more practical 
for molecular dynamics simulations especially for finite 
temperatures. Moreover, when analytical potentials are 
available, the effect of various types of interactions can be 
investigated by switching off, or modifying, correspond- 
ing potential energy terms, being able to improve in this 
way our understanding of the studied problem [18| . The 
obtained empirical potentials are then used to calculate, 
through molecular dynamics simulations, the mechani- 
cal response of graphene under various loads. To test 
these results, density functional theory calculations of 
graphene's elastic behavior have been also performed. 

The derived bond stretching and angle bending force 
fields are discussed in the Section [Til along with the an- 
alytical formulas for the respective potential energies. 
Then, the Section IIIII contains results concerning the 
mechanical response of graphene monolayers on uniax- 
ial tension, hydrostatic tension, and shear stress. Stress- 
strain curves and the corresponding elastic parameters 
are presented. Finally, the Section |IV] concludes our 
work. 



II. IN-PLANE FORCE FIELDS 

We performed first principles', Density Functional 
Theory (DFT) calculations at the GGA/PBE levelQl] 
for an infinite graphene sheet. We used the Quantum- 




FIG. 1: Bond stretching potential energy as a function of 
the interatomic distance between bonded carbon atoms in 
graphene layers. Circles show the results obtained through 
density functional theory. Solid line represents the best fit 
with a Morse potential, Eq. fTJ). 



FIG. 2: Bending energy per hexagonal ring as a function of 
the angle deformation shown in the inset. Circles show the 
results obtained through density functional theory. Solid line 
represents the best fit with Eq. <(3j , from where the parame- 
ters of the angle bending potential ((2]) are derived. 



Espresso computer code[20| with a ultra-soft pseudopo- 
tential [2l[ generated by a modified RRKJ approach [22[ 
which is proven to reproduce accurately structural, vi- 
brational and thermodynamic properties of carbon al- 
lotropes (23|. Finally, we adopted a minimal two-atom 
unit cell, a 10x10 k-mesh and cutoffs 70 and 560 Ry for 
the wave functions and charge density respectively. 

To find the energy dependence of the bond stretch- 
ing motion we scaled the lattice, thus stretching all C- 
C bonds without altering any valence angle. In Fig. Q] 
the calculated energy per bond is shown. For efficient 
atomistic simulations, we have fitted the obtained DFT 
results with appropriate anharmonic functions, thus pro- 
viding analytical potentials for the description of the in- 
plane molecular mechanics of graphene. Regarding bond 
stretching, the Morse potential 
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has been considered. Through fitting of the numerical 
DFT results with eq. (1) , the parameter values D = 
5.7 eV, a = 1.96 A^ 1 , rg = 1.42 A have been obtained 
(see Fig. [I]). 

The angle bending force field is described by a nonlin- 
ear potential containing quadratic and cubic terms 
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To obtain the parameters k and k! of this expression, we 
consider the deformed graphene lattice where each ben- 
zene ring is identically deformed as shown in the inset 
of Fig. [2] only the valence angles have been uniformly 



altered by x or x/2, while all bond lengths have their 
equilibrium value of 1.42 A. In the same figure we show 
the dependence of the DFT calculated total energy per 
ring on the angle-bending parameter x. For such a defor- 
mation, considering the angle bending potential energy 
of Eq. @, the obtained energy per ring is 
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(3) 

Fitting the DFT results of Fig.[2]with the last expression 
yields the parameter values k = 7.0 eV/rad 2 and k' = 
4 eV/rad 3 for the angle bending potential, Eq.([2]). 

We mention that bond stretching and angle bending 
potentials like those discussed here, Eqs. ((TJ) and ©, 
have been also presented in Ref. [24|, but with differ- 
ent parameter values: D = 234.42 kcal/mol = 10.1 eV, 
a = 1.75 A" 1 , r = 1.27 A, k = 36.26 kcal/mol/rad 2 = 
1.56 eV/rad 2 , and k' = (a linear angle bending poten- 
tial has been considered in that work). 



III. MECHANICAL RESPONSE OF 
GRAPHENE MONOLAYERS 

Stress-strain curves and elastic constants of graphene 
monolayers have been calculated using both first prin- 
ciple's method and molecular dynamics (MD) atomistic 
simulations with the potentials presented in the previous 
section. 

In the former case, DFT calculations were performed, 
as previously, but for a larger rectangular unit cell com- 
prising by 16 atoms and a 6x12 k-mesh. Strains were 



applied in one (either the zigzag or armchair, for uniax- 
ial tensions) or both (for hydrostatic tension) of the two 
vertical directions of the unit cell. For a given strain level 
the corresponding unit cell dimensions are fixed while all 
other structural parameters such as atomic positions and 
vertical to strain unit-cell dimension for uniaxial strain, 
were allowed to relax. The stresses (forces per length) 
were calculated as numerical derivatives of the total en- 
ergy with respect to the appropriate length (or area) of 
the unit cell for given strain. 

In order to obtain the mechanical response of graphene 
at a fixed force, /, acting on one or more of its edges, 
atomistic simulations have been performed as follows: we 
start with the equilibrium structure of graphene, with- 
out any external force (i.e. all first neighboring distances 
are 1.42 A and all valence angles 2tt/3). Then, a con- 
stant force / is applied at all atoms of the appropriate 
edge, depending on the case (the force may be perpen- 
dicular or parallel to the edge). Now, the system is out 
of equilibrium and Newton's equations of motion are nu- 
merically solved, applying a friction term at each atom, 
to follow the evolution of the system. Then, the sys- 
tem gradually goes to a new equilibrium, compatible to 
the applied forces. In this process, due to the friction 
the total energy of the system decreases from its initial 
value towards a new lower value, corresponding to the 
new equilibrium (the deformed state due to the applied 
stress). The total kinetic energy increases initially from 
zero and after some decaying oscillations it gradually van- 
ishes, when the new equilibrium has been reached. Dur- 
ing this evolution, strains are developed in the graphene 
sheet and, following a transient oscillatory behavior, they 
converge to the equilibrium values corresponding to the 
applied stress. To avoid finite size effects, the equilibrium 
strains at the middle of the examined graphene sheets are 
recorded. Then the same procedure is repeated for an- 
other value of the applied force /. Graphene monolayers 
consisting of 7482 and 17030 carbon atoms have been 
simulated, providing the same results. The friction co- 
efficient was 10 psec -1 . It has been checked that other 
values of the friction coefficient result in the same de- 
formed state, affecting only the transient period and the 
number of oscillations needed before reaching the corre- 
sponding equilibrium. 



A. Uniaxial tension 

Here, atomistic simulations have been performed, as 
described above, with a constant force / applied at all 
atoms belonging to the two opposite edges of graphene. 
The force is perpendicular to the corresponding edges and 
it is directed outwards. Fixed strains at the two opposite 
edges of the unit cell have been imposed in the respective 
DFT calculations, while the vertical direction has been 
relaxed to minimize the total energy. For a given strain, 
the stress is obtained through the numerical derivative 
of the total energy with respect to the size of the unit 
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FIG. 3: Stress-strain curves for uniaxial deformations of 
graphene monolayers. Top: The uniaxial tension is applied 
to graphene's zigzag edges. Bottom: The uniaxial tension is 
applied to graphene's armchair edges. Filled circles show re- 
sults obtained through molecular dynamics simulations, while 
open squares correspond to density functional theory calcu- 
lations. Solid lines represent fits with the nonlinear relation, 
Eq. 0. The insets depict the variation of the Poisson ratio 
with the uniaxial strain. 



cell along the direction where the strains are imposed. 
MD (DFT) results are presented for the cases where the 
tensile force (strain) is applied either on the zigzag or on 
the armchair edges of graphene. 

Fig.Odepicts the force per unit length (the correspond- 
ing stress at the edge of a two-dimensional material) as a 
function of the strain (%) for uniaxial tensions applied on 
the graphene's zigzag (top) or armchair (bottom) edges. 
In the former case, the DFT and MD results coincide up 
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to a strain of 14-15%, above of which failure is obtained 
in the atomistic simulations. In the latter case, the DFT 
and MD calculations are almost identical up to about 
12-14% strain, while failure is obtained at a strain larger 
than 28% in the atomistic simulations and a bit earlier 
in the first principles method. 

The slope of the obtained stress-strain curves in the 
limit of small strains and stresses leads to a 2D Young 
modulus E2D = 320 N/m, for uniaxial tensions applied 
on both types of graphene edges. Considering the thick- 
ness of a graphene monolayer to be I — 0.335 nm (the 
interlayer spacing of graphite), the derived effective 3D 
Young modulus is E eff = E 2D /l = 0.96 TPa. The 2D 
intrinsic strengths (the maximum tensile stresses sup- 
ported before failure) are a 2D = 32 — 34 N/m for uniaxial 
forces applied on the zigzag edges and a 2D = 39—45 N/ m 
for uniaxial forces applied on the armchair edges. The ef- 
fective 3D intrinsic strengths, <J e ff — &2d/1, are around 
100 GPa and 120-130 GPa, respectively. The failure 
observed in atomistic MD simulations for stresses larger 
than the intrinsic strength, occurs as shown in Fig. 7 of 
Ref. [25^ . i.e. when the tensile force is applied on the 
armchair (zigzag) edges of graphene, the terminal single 
layer of atoms on which the forces are applied (the ter- 
minal layer of hexagonal atomic rings) breaks and it is 
disrupted from the system. 

The calculated Young modulus and intrinsic strengths 
are in agreement with available experimental estimates 
and in the range of theoretically computed values from 
other works. Early measurements of the elastic constants 
of highly ordered pyrolytic graphite indicated a value 
of 1 TPa for the in-plane Young modulus. More re- 
cently, indirect experimental estimates, obtained through 
nanoindentation measurements of free-standing mono- 
layer graphene result in 2D Young modulus E 2 d = 
340 ± 50 N/m and strength a 2D = 42 ± 4 N/m, lead- 
ing to effective 3D values E e ff = 1.0 ± 0.1 TPa and 
a e ff = 130 ± 20 GPa, respectively. Theoretical pre- 
dictions for the Young modulus of graphene have been 
obtained through a variety of methods. Molecular me- 
chanics simulations yield E eff = 0.95 TPa [Hj], MD 
simulations E e ff — 1.01 TPa @, Monte Carlo cal- 
culations E ef} = 1.04 TPa (E 2D = 350 N/m) [|, 
tight-binding approximations E e ff = 0.91 TPa 0, and 
density functional theory methods result in values of 
E eff = 1.05 - 1.09 TPa SUS^. Intrinsic strengths 
<j z e ^ = 110 GPa and a^jf = 121 GPa for forces applied 
along the zigzag and armchair edges, respectively, have 
been calculated in Ref. [6j. 

The elastic response of monolayer graphene under uni- 
axial extension can be captured by the following nonlin- 
ear relation 

C2D = E 2D ■ e + D 2D ■ e 2 , (4) 

where e is the uniaxial strain and D 2 o is the third- 
order elastic modulus which is typically negative [l3| . 
The Young modulus is about 320 N/m for both direc- 
tions of uniaxial tension and the D 2 £> is derived through 



fitting of the numerical data with Eq. When ten- 
sion is applied on the zigzag edges, the extracted D 2 jj is 
-700 N/m in the DFT and -670 N/m in the MD re- 
sults. For the armchair direction the obtained D 2 d values 
are —670 N/m and —560 N/m, respectively. These val- 
ues are in agreement with the experimentally determined 
value of-690 iV/m [H[. 

The insets in Fig. [3] show the corresponding Poisson's 
ratio, v , as a function of the uniaxial strain. In both di- 
rections, a value of v ss 0.22 is obtained for small strains, 
while the Poisson ratio decreases for larger strains. When 
the uniaxial tension is applied to the zigzag edges, the 
slope of the almost linear decrease of the Poisson ra- 
tio is -0. 006 /%strain in both MD and DFT calcula- 
tions. In the other direction, atomistic simulations give 
a linear decrease of v with a slope of — 0.007/%siram, 
while the DFT calculations coincide with the MD ones 
for strains up to about 12%, but the slope changes for 
larger strains resulting in an approximate average slope 
of — 0.004/%siram. Similar behavior has been obtained 
from the DFT results of Ref. @, where v starts from 
a value around 0.19 at small strains and then it al- 
most linearly decreases up to strains 25-30%, exhibiting 
a more abrupt slope in the case of tension applied on the 
zigzag edges. Other theoretical estimates give values of 
v around 0.15 at °K H, v = 0.28 at 1 atm and 300 °K 
[13, and v = 0.21 at 300 °K [|. 



B. Hydrostatic tension 

In this case, constant forces are applied at all boundary 
atoms of the graphene sheets studied with MD, in such a 
way that the corresponding stress (force per unit length) 
is the same at both kinds of edges. Taking into account 
that the average interatomic distances along the zigzag 
and armchair edges are V3ro/2 and 0.75ro, respectively, 
the applied forces have different magnitudes at each type 
of edge. They are perpendicular to the edges and are 
directed outwards. In the DFT calculations, a uniform 
strain in both vertical directions was applied. The hydro- 
static stress was calculated by numerical differentiation 
of the energy with respect to the unit cell's area, at given 
uniform strains. Applying uniform strain in both direc- 
tions, i.e. assuming that the material is isotropic, is a 
reasonable approximation as it is justified below. 

Fig [4] presents the force per unit length as a function of 
the relative surface change DS/Sq = (S — Sq)/Sq, where 
So is the initial undeformed surface of the system and S 
is the final equilibrium surface corresponding to the ap- 
plied forces. Graphene supports hydrostatic tension up 
to more than 30% relative change of its surface. The max- 
imum tensile hydrostatic stress before failure is obtained 
around 31-32 N/m, which, divided by I, corresponds to 
an effective 3D stress 90-100 GPa. This value of maxi- 
mum stress is close to the intrinsic strength of uniaxial 
tension applied on the zigzag edge, which is the minimum 
of the two values of intrinsic strengths corresponding to 
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FIG. 4: Stress-strain curve for hydrostatic tension defor- 
mation of graphene monolayers. Filled circles (open squares) 
show results obtained through molecular dynamics simula- 
tions (density functional theory calculations). The inset de- 
picts the force per unit length as a function of the relative 
length deformation along the armchair (filled circles) and the 
zigzag (open squares) direction, obtained through the molec- 
ular dynamics calculations. Dashed lines are guides to the 
eye. 



the two different directions of uniaxiai tension. How- 
ever, this does not imply that hydrostatic failure occurs 
at the zigzag edges by the same mechanism as in the cor- 
responding uniaxial tension 29]. On the contrary, MD 
simulations indicate a more complicated failure pattern 
in this case. 

In the small strain/stress limit, the slope of the stress- 
strain curve shown in Fig |4] results in a 2D bulk modulus 
B 2D = 200 N/m. This value of the 2D bulk modulus 
is in agreement with estimates of B^d ~ 13 eV/A 2 « 
200 N/m, obtained from Monte Carlo simulations |8j. 

In order to investigate potential asymmetries in the 
elastic response of graphene under hydrostatic load, we 
show in the inset of Fig|4]the relation between the applied 
force per unit length and the relative extensions along the 
zigzag and armchair directions, separately, as obtained 
from the atomistic simulations. A completely symmetric 
response arises, as the strains along the two directions 
are identical. This confirms the equivalence between the 
DFT calculations, obtained by applying uniform strains 
along the two directions, and the hydrostatic MD results, 
where uniform stresses are applied in the two directions. 
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FIG. 5: Stress-strain curves for shear deformations of 
graphene monolayers, obtained through molecular dynamics 
simulations. Filled circles (open squares) show results corre- 
sponding to shear forces acting on the armchair edge (zigzag 
edge) of graphene. Dashed lines are guides to the eye. 



ing the whole MD process, while constant forces are ap- 
plied at all atoms of the opposite edge. The direction 
of the forces is parallel to the edge on which they ap- 
ply. The shear response of both kinds of edges has been 
calculated through this procedure. 

Fig. [5] shows the shear stress as a function of shear 
strain (the tangent of the strain angle at the middle of 
the examined graphene sheet). A noticeable asymmetry 
on the shear response of the zigzag and armchair edges 
is obtained. The armchair edges can withstand shear 
stresses up to around 11 N/m, where a maximum shear 
strain 0.14 is developed. A significantly smaller maxi- 
mum shear stress of less than 4 N/m is supported by the 
zigzag edges, but with a larger shear strain. 

Almost linear stress-strain curves describe the shear 
deformations of both zigzag and armchair edges of 
graphene. The obtained 2D shear moduli are G% D = 
84 N/m for shear stresses applied on the armchair edge 



and G Z 2D = 22 N/m for shear applied on the zigzag edge. 
The effective 3D shear moduli G e // = G2D/I are 250 
GPa and 65 GPa, respectively. The asymmetric shear 
response of the two kinds of edges result in a difference 
by a factor of about 4 of the corresponding elastic mod- 
uli. A higher value of shear modulus in graphene, around 
150 N/m, has been obtained from Monte Carlo simula- 
tions Hi. 



C. Shear stress 

Shear deformations are investigated merely by MD 
simulations. To this end, one edge of graphene has been 
kept fixed, at its free-of-load equilibrium structure, dur- 



D. Effect of the anharmonicity of force fields 

As it has been mentioned in the Sectionfl] one of the ad- 
vantages of atomistic simulations, when analytical force 
fields are available, is that the influence of several fac- 
tors can be explored by altering or eliminating respective 



potential energy terms. Here, we examine the effect of 
nonlinear interactions on the uniaxial tensile response of 
graphene, since both bond stretching and angle bending 
potentials, Eqs. ([1]) and @, presented in this work are 
anharmonic. 

Using respective MD simulations, results have been 
obtained for the cases that (i) the angle bending po- 
tential is linearized by ignoring the cubic term (setting 
k' = 0) in Eq. ([2]) and (ii) both angle bending and bond 
stretching potentials are linearized; i.e. additionally to 
the change mentioned previously, the Morse potential is 
also linearized through the substitution of Eq. ([TJ) by its 
harmonic approximation V* mear — Da 2 (r — ro) 2 . 

Fig. [B] compares the stress-strain curves on uniaxial 
tension obtained in the cases (i) and (ii) with the corre- 
sponding results shown in Fig. [31 where the anharmonic 
potentials of Eqs. ([T]) and ([2]) have been used. The in- 
sets present these comparisons regarding the strain de- 
pendence of the Poisson ratio. It can be seen that the 
nonlincarity of the angle bending potential has not any 
noticeable effect on the calculated mechanical response. 
Instead, the response is drastically different when har- 
monic bond stretching potentials are considered too. As 
expected, in the latter case a linear stress-strain relation 
is obtained (deviating from the nonlinear response for 
strains above a few per cent) and no failure is present. 
Further, the Poisson ratio shows no dependence on the 
strain, in marked contrast to the behavior observed for 
nonlinear bond stretching force fields. 



IV. CONCLUSIONS 

Empirical bond stretching and angle bending poten- 
tials have been presented, suitable for atomistic simula- 
tions of graphene. The potentials have been derived using 
methods from first principles and are expected to accu- 
rately represent the respective interactions in graphene. 
Appropriate nonlinear functions that analytically de- 
scribe the numerically obtained potentials are provided, 
for efficient use in molecular mechanics calculations. 

The mechanical response of graphene under various 
loads is examined, through both molecular dynamics 
simulations, using the calculated empirical force fields, 
and density functional theory calculations. Stress-strain 
curves have been presented for uniaxial tensions or shear 
stresses applied on both armchair and zigzag edges of 
graphene, as well as for hydrostatic tension. The de- 
pendence of Poisson ratio on strain is also calculated in 
the case of uniaxial tensions. Corresponding elastic mod- 
uli, such as Young modulus, two-dimensional bulk mod- 
ulus, or shear moduli, and maximum values of supported 
stresses or strain before failure are discussed and com- 
pared with experimental estimates, when available, and 
other theoretical results. 

Finally, the effect of the nonlinearity of the potentials 
on the uniaxial tensile response of graphene is consid- 
ered. We find that the anharmonicity in the angle bend- 
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FIG. 6: The effect of the anharmonicity of the potentials on 
the uniaxial deformations of graphene monolayers. Filled cir- 
cles show results using the fully nonlinear potentials of Eqs. 
([T]) and open squares correspond to linearized angle bend- 
ing potential and nonlinear bond stretching potential, while 
open diamonds correspond to completely linearized poten- 
tials. All the results are obtained using molecular dynam- 
ics simulations. Top: The uniaxial deformation is applied to 
graphene's zigzag edges. Bottom: The uniaxial deformation 
is applied to graphene's armchair edges. The insets depict 
the corresponding results for the Poisson ratio as a function 
of the uniaxial strain. Dashed lines are guides to the eye. 



ing potential plays no role in the examined mechanical re- 
sponse, but, instead, nonlinear terms of the bond stretch- 
ing potential crucially affect the obtained elastic behavior 
and, thus, should be necessarily taken into account. 
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